home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Digital / ZSolve.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  10.0 KB  |  311 lines
  1. (*  :Title:    Solution of Difference Equations via z-transforms. *)
  2.  
  3. (*  :Authors:    Wally McClure, Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    This package will solve linear constant coefficient
  7.         difference equations by means of the z-transform.
  8.  *)
  9.  
  10. (*  :Context:    SignalProcessing`Digital`ZSolve`  *)
  11.  
  12. (*  :PackageVersion:  2.7    *)
  13.  
  14. (*
  15.     :Copyright:    Copyright 1990-1991 by Brian L. Evans
  16.         Georgia Tech Research Corporation
  17.  
  18.     Permission to use, copy, modify, and distribute this software
  19.     and its documentation for any purpose and without fee is
  20.     hereby granted, provided that the above copyright notice
  21.     appear in all copies and that both that copyright notice and
  22.     this permission notice appear in supporting documentation,
  23.     and that the name of the Georgia Tech Research Corporation,
  24.     Georgia Tech, or Georgia Institute of Technology not be used
  25.     in advertising or publicity pertaining to distribution of the
  26.     software without specific, written prior permission.  Georgia
  27.     Tech makes no representations about the suitability of this
  28.     software for any purpose.  It is provided "as is" without
  29.     express or implied warranty.
  30.  *)
  31.  
  32. (*
  33.     :History:    Date started -- April 30, 1990
  34.         Integrated into signal processing packages -- June 6, 1990
  35.         Handle arbitrary initial conditions -- May 9, 1991
  36.  
  37.         This package started up from a special topics class.  
  38.  *)
  39.  
  40. (*  :Keywords:    linear constant-coefficient difference equations    *)
  41.  
  42. (*
  43.     :Source:    {Fundamentals of Digital Signal Processing}
  44.           by Lonnie C. Ludeman 
  45.  *)
  46.  
  47. (*  :Warning:    *)
  48.  
  49. (*  :Mathematica Version:  1.2 or 2.0  *)
  50.  
  51. (*  :Limitation:  *)
  52.  
  53. (*  :Discussion:  *)
  54.  
  55. (*  :Functions:   PartialZ  ZSolve  *)
  56.  
  57.  
  58.  
  59. (*  B E G I N     P A C K A G E  *)
  60.  
  61. BeginPackage[ "SignalProcessing`Digital`ZSolve`",
  62.           "SignalProcessing`Digital`ZTransform`",
  63.           "SignalProcessing`Digital`InvZTransform`",
  64.           "SignalProcessing`Digital`ZSupport`",
  65.           "SignalProcessing`Support`TransSupport`",
  66.           "SignalProcessing`Support`ROC`",
  67.           "SignalProcessing`Support`SigProc`",
  68.           "SignalProcessing`Support`SupCode`" ]
  69.  
  70.  
  71. If [ TrueQ[ $VersionNumber >= 2.0 ],
  72.      Off[ General::spell ];
  73.      Off[ General::spell1 ] ];
  74.  
  75.  
  76. (*  U S A G E     I N F O R M A T I O N  *)
  77.  
  78. RightSided::usage =
  79.     "RightSided is an option for ZSolve. \
  80.     If True, then the solution to the difference equation will be \
  81.     right-sided, left-sided otherwise."
  82.  
  83. ZSolve::usage =
  84.     "ZSolve[ diffequ == drivingfun, y[n] ] solves the difference \
  85.     equation diffequ = drivingfun, where diffequ is a linear constant \
  86.     coefficient difference equation and drivingfun is the driving \
  87.     function (a function of n). \
  88.     Thus, diffequ has the form a0 y[n] + a1 y[n - 1] + .... \
  89.     By using options, one can specify initial values; e.g., \
  90.     ZSolve[ y[n] + 3/2 y[n-1] + 1/2 y[n-2] == (1/4)^n Step[n], \
  91.     y[n], y[-1] -> 4, y[-2] -> 10 ]. \
  92.     A difference equations of N terms needs N-1 initial values. \
  93.     All unspecified initial values are considered to be zero. \
  94.     Right-sided and left-sided solutions are possible. \
  95.     ZSolve can justify its answers."
  96.  
  97. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  98.  
  99.  
  100. Begin[ "`Private`" ]
  101.  
  102.  
  103. (*  Messages  *)
  104. ZSolve::notsupported =
  105.     "You specified initial conditions that start too far back."
  106. ZSolve::noztransform =
  107.     "The forcing function `` does not have a valid z-transform."
  108. ZSolve::redundant =
  109.     "Same initial index is assigned two different values."
  110. ZSolve::toomany =
  111.     "`` initial conditions given for a difference equation of order ``."
  112. PartialZ::noztransform =
  113.     "The difference equation `` does not have a partial z-transform."
  114.  
  115.  
  116. (* introduction *)
  117. introduction [ y_[n_], minindex_, maxindex_, diffequ_,
  118.         minoffset_, maxoffset_, drivingfun_, initlist_ ] :=
  119. Block [    {rstr},
  120.     Print["Solving for ", y[n], " in the difference equation"];
  121.     Print[ "  ", diffequ, " = ", drivingfun ];
  122.     rstr = If [ TrueQ[rightsided],
  123.             StringForm["`` > ``", n, maxindex],
  124.             StringForm["`` < ``", n, minindex] ];
  125.     Print[ "  such that ", rstr,
  126.            " and subject to the initial conditions" ];
  127.     If [ EmptyQ[initlist],
  128.          Print[ "  being zero." ];
  129.          Print[ "  ", initlist ] ];
  130.     Print[ " " ];
  131.     Print[ "The first step is to substitute the initial conditions" ];
  132.     Print[ "into the difference equation for ", n, " from ",
  133.         minindex - maxoffset, " to ",
  134.         maxindex - maxoffset, "." ];
  135.     Print[ "This will give rise to ", maxindex - minindex + 1,
  136.         " impulse function(s) which will" ];
  137.     Print[ "be added to the right-hand side of the difference equation." ];
  138.     Print[ "Then, the difference equation will be valid for all values" ];
  139.     Print[ "of ", n, ", so the unilateral z-transform can used to solve" ];
  140.     Print[ "the problem." ];
  141.     Print[ " " ] ]
  142.  
  143.  
  144. (*  homogeneous  *)
  145. homogeneous[ drivingfun_, diffequ_, y_[n_], dialogue_, rsided_, {} ] :=
  146.     diffequ == drivingfun
  147.  
  148. homogeneous[ drivingfun_, diffequ_, y_[n_], dialogue_, rsided_, initlist_ ] :=
  149.     Block [    {auginitlist, homtable, indexlist, initcond, matrix, maxindex,
  150.          minindex, maxoffset, minoffset, newdiffequ, n0,
  151.          ntemp, offsetlist, order, R, step},
  152.  
  153.         (* determine the order of the difference equation *)
  154.         offsetlist = Apply[List, diffequ] /.
  155.                  {a_ y[n + b_] :> b, y[n + c_] :> c, y[n] -> 0};
  156.         offsetlist = Select[offsetlist, IntegerQ];
  157.         maxoffset = Max[offsetlist];
  158.         minoffset = Min[offsetlist];
  159.         order = maxoffset - minoffset;
  160.  
  161.         (* determine the maximum index of the initial conditions *)
  162.         If [ EmptyQ[initlist],
  163.              minindex = maxindex = 0,
  164.  
  165.              indexlist = ReplaceAll[ initlist,
  166.                          (y[index_] -> value_) :> index ];
  167.              If [ ! SameQ[indexlist, Union[indexlist]],
  168.               MyMessage[ZSolve::redundant, Null] ];
  169.              initcond = Length[indexlist];
  170.              If [ initcond > order,
  171.               Message[ZSolve::toomany, initcond, order];
  172.               indexlist = Take[indexlist, -order];
  173.               initcond = order ];
  174.              maxindex = Max[indexlist];
  175.              minindex = Min[indexlist];
  176.              If [ ( maxindex - minindex ) >= order,
  177.               Message[ZSolve::notsupported] ] ];
  178.  
  179.         (* update dialogue with user *)
  180.         If [ SameQ[dialogue, All],
  181.              introduction[ y[n], minindex, maxindex, diffequ,
  182.                 minoffset, maxoffset, drivingfun, initlist ] ];
  183.  
  184.         (* determine effects of initial conditions *)
  185.         newdiffequ = diffequ /. n -> ntemp;
  186.         homtable = Table[(newdiffequ) Impulse[n - ntemp],
  187.                  {ntemp, minindex - maxoffset,
  188.                      maxindex - maxoffset}];
  189.         auginitlist = Append[initlist, y[index_] :> 0];
  190.         step = If [ TrueQ[rightsided],
  191.                 Step[n - (maxindex + 1)],
  192.                 Step[(mindex - 1) - n ] ];
  193.  
  194.         diffequ == drivingfun step +
  195.                 Apply[Plus, homtable /. auginitlist] ]
  196.  
  197.  
  198. (*  preprocess  *)
  199. preprocess[ left_ == right_, y_[n_], dialogue_, rs_, initlist_ ] :=
  200.     Block [    {diffequ, drivingfun},
  201.  
  202.         diffequ = Select[left - right,
  203.                  MatchQ[#1, a_. y[n + n0_.]]&];
  204.         drivingfun = - Select[left - right,
  205.                       ! MatchQ[#1, a_. y[n + n0_.]]&];
  206.  
  207.         homogeneous[drivingfun, diffequ, y[n], dialogue, rs, initlist] ]
  208.  
  209.  
  210. (*  PartialZ  *)
  211. PartialZ[ x_ + v_, rest__ ] :=
  212.     PartialZ[ x, rest ] + PartialZ[ v, rest ]
  213.  
  214. PartialZ[ a_ x_, y_[n_], z_, YY_ ] :=
  215.     a PartialZ[ x, y[n], z, YY ] /; FreeQ[a, n]
  216.  
  217. PartialZ[ y_[n_ + n0_.], y_[n_], z_, YY_ ] :=
  218.     z^n0 YY /; FreeQ[n0, n]
  219.  
  220. PartialZ[ a_, y_[n_], z_, YY_ ] :=
  221.     MyMessage[ PartialZ::noztransform, 0, a ]
  222.  
  223.  
  224. (*  ZSolve  *)
  225. ZSolve/: Options[ZSolve] := { Dialogue -> True, RightSided -> True }
  226.  
  227. ZSolve[ y_[n_ + a_.] == drivingfun_, y_[n_] ] :=
  228.     ( drivingfun /. { n -> n - a } )
  229.  
  230. ZSolve[ left_ == right_, y_[n_], init___ ] :=
  231.     Block [    {dialogue, diffequ, drivingfun, equationlist, index,
  232.          initlist, options, rightsided, rminus, rplus, X, Y,
  233.          Ycollected, Ypartial, YY, YYstr, z, zobj, zsol},
  234.  
  235.         (* Parse options, determine the level of Dialogue requested *)
  236.         options = ToList[init] ~Join~ Options[ZSolve];
  237.         dialogue = Replace[Dialogue, options];
  238.         rightsided = Replace[RightSided, options];
  239.         initlist = Select[ options,
  240.                    MatchQ[#1, y[index_] -> value_]& ];
  241.  
  242.         (* Collect terms, augment driving function with i.c.'s *)
  243.         equationlist = preprocess[ left == right, y[n],
  244.                        dialogue, rightsided, initlist ];
  245.         diffequ = First[equationlist];
  246.         drivingfun = Second[equationlist];
  247.  
  248.         (* Tell the user about the augmented difference equation *)
  249.         If [ dialogue ||  SameQ[dialogue, All],
  250.              Print["Solving for ", y[n], " in the difference equation"];
  251.              Print["augmented by the initial conditions:"];
  252.              Print["  ", diffequ, " = ", drivingfun] ];
  253.  
  254.         (* Find complete z-transform of driving function *)
  255.         zobj = ZTransform[ drivingfun, n, z, Dialogue -> False ];
  256.         If [ ! SameQ[Head[zobj], ZTransData],
  257.              Message[ ZSolve::noztransform, drivingfun ];
  258.              Return[] ];
  259.         rminus = GetRMinus[zobj];             (* R- of X(z) *)
  260.         rplus = GetRPlus[zobj];                 (* R+ of X(z) *)
  261.         X = TheFunction[zobj];             (* X(z) from z object *)
  262.  
  263.         (* Find partial z-transform of difference equation *)
  264.         Ypartial = PartialZ[ diffequ, y[n], z, YY ];
  265.  
  266.         (* Dialogue taking z-transform of both sides to user *)
  267.         If [ dialogue,
  268.              Print[ "After taking the z-transform of both sides and" ];
  269.              Print[ "solving for the z-transform of ", y[n], ":" ] ];
  270.         If [ SameQ[dialogue, All],
  271.              Ycollected = Collect[Ypartial, YY];
  272.              YYstr = ToString[StringForm["Z{ `` }", y[n]]];
  273.              Print[ "The z-transform of the left side is:" ];
  274.              Print[ "  ",
  275.                     ( Ycollected /. YY -> YYstr ) /. z -> "z" ];
  276.              Print[ "The z-transform of the right side is:" ];
  277.              Print[ "  ", X /. z -> "z" ];
  278.              Print[ "Solving for the unknown z-transform yields" ] ];
  279.  
  280.         (* Inverse transform complete z-transform of y, Y(z) *)
  281.         zsol = Solve[ ( Ypartial == X ) /. z -> z^-1, YY ];
  282.         Y = ( YY /. First[ zsol ] ) /. z -> z^-1;
  283.  
  284.         If [ dialogue || SameQ[dialogue, All],
  285.              Print[ "  ", Y /. z -> "z" ];
  286.              Print[ "Inverse transforming this gives ", y[n], ":" ] ];
  287.  
  288.         InvZTransform[{Y, rminus, rplus}, z, n, Dialogue -> False] ]
  289.  
  290.  
  291. (*  E N D     P A C K A G E  *)
  292.  
  293. End[]
  294. EndPackage[]
  295.  
  296. If [ TrueQ[ $VersionNumber >= 2.0 ],
  297.      On[ General::spell ];
  298.      On[ General::spell1 ] ];
  299.  
  300.  
  301. (*  H E L P     I N F O R M A T I O N  *)
  302.  
  303. Combine[ SPfunctions, { ZSolve } ]
  304. Protect[ ZSolve ]
  305.  
  306.  
  307. (*  E N D I N G     M E S S A G E  *)
  308.  
  309. Print[ "ZSolve, a difference equation solver, is loaded." ]
  310. Null
  311.